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ON . We consider tunneling of electromagnetic waves through a polariton band gap of a 1-D chain of 

ON ■ atoms. We analytically show that a defect embedded in the structure gives rise to the resonance 

transmission at the frequency of a local polariton state associated with the defect. Numerical 
Monte-Carlo simulations are used to examine properties of the electromagnetic band arising inside 
the polariton gap due to finite concentration of defects. 
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£ . I. INTRODUCTION 

The resonant tunneling of electromagnetic waves through different types o£,|Optical barriers is a fast developing 
area of optical physics. This effect was first considered for photonic crystals, lira where forbidden band-gaps in the 
electromagnetic speeiaim form optical barriers. Macroscopic defects embedded in the photonic crystal give rise to 
local photon modes|3i3 which induce the resonant transmission of electromagnetic waves through the band-gaps. 

A different type of photonic band gaps arises in polar dielectrics, where a strong resonance interaction between 
the electromagnetic field and dipole active internal excitations of a dielectric brings about a gap between different 
branches of polaritons. Recently pit-was suggested that regular microscopic impurities embedded in such a dielectric 
give rise to local polariton states,ETt2l where a photon is_cpupled to an intrinsic excitation of a crystal, and both these 
components are localized in the vicinity of the defect. lid The main peculiarity of the local polaritons is that their 
electromagnetic component is bound by a microscopic defect whose size is many order of magnitude smaller then the 
wavelengths of respective photons. Another important property of these states is the absence of a threshold for their 
appearance even in 3-D isotropic systems, while for all other known local states the "strength" of a defect must exceed 
a certain critical value before the state would split off a continous spectrum. The reason for this peculiar behavior is 
a strong van Hove singularity in the polariton density of states in the long wave region, caused by a negative effective 
mass of the polariton-forming excitations of a crystal. 

The feasibility of resonant electromagnetic tunneling induced by local polaritons, however, is not self-evident. The 
idea of a microscopic defect affecting propagation of waves with macroscopic wavelength seems to be in contradiction 
ON ' with common wisdom. Besides, it was found that the energy of the electromagnetic component of local polaritons 
is very small compared to the energy of its crystal counterpart. The existence of the tunneling effect was first 
numerically demonstrated in Ref. ^, where a 1-D chain of dipoles interacting with a scalar field imitating transverse 
electromagnetic waves was considered. It was found that a single defect embedded in such a chain results in near 
100% transmission at the frequency of local polaritons through a relatively short chain (50 atoms). The frequency 
profile of the transmission was found to be strongly assymetric, in contrast to the case of electron tunneling£j 

In most cases (at least for a small concentration of the transmitting centers) one-dimensional models give a reliable 
description of tunneling processes, because by virtue of tunneling, a wave ^propagates along a chain of resonance 
O ■ centers, for which a 1-D topology has the highest probability of occurrence.tll In our situation, it is also important 
that the local polariton states (transmitting centers) occur without a threshold in 3-D systems as well as in 1-D 
systems. This ensures that the transmission resonances found in Ref. ^ are not artifacts of the one-dimensional nature 
of the model, and justifies a further development of the model. In the present paper we pursue this development in 
two interconnected directions. First, we present an exact analytical solution of the transmission problem through the 
chain with a single defect. This solution explains the unusual asymmetric shape of the transmission profile found 
in numerical calculationsLi and provides insight into the phenomenon under consideration. Second, we carry out 
numerical Monte-Carlo simulation of the electromagnetic transmission through a macroscopically long chain with a 
finite concentration of defects, and study the development of a defect-induced electromagnetic pass band within the 
polariton band-gap. The analytical solution of a single-defect model allows us to suggest a physical interpretation for 
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some of the peculiarities of the transmission found in numerical simulations. As a by-product of our numerical results 
we present a new algorithm used for the computation of the transmission--. This algorithms is based upon a blend 
of transfer-matrix approach with ideas of the invariant-embedding method,E3 and proves to be extremely stable even 
deep inside the band-gap, where traditional methods would not work. 

Though we consider the one-dimensional model, the results obtained are suggestive for experimental observation of 
the predicted effects. Actually the damping of the electromagnetic waves is more experimentally restrictive than the 
topology of the system. We, however, discuss the effects due to damping and come to the conclusion that the effects 
under discussion can be observed in regular ionic crystals in the region of their phonon-polariton band-gaps. 

The paper is organized as follows. The Introduction is followed by an analytical solution of the transmission problem 
in a single-impurity situation. The next section presents results of Monte-Carlo computer simulations. The algorithm 
used in numerical calculations is derived and discussed in the Appendix. The paper concludes with a discussion of 
the results. 



II. DESCRIPTION OF THE MODEL AND ANALYTICAL SOLUTION OF A SINGLE-DEFECT 

PROBLEM 



1. The model 



Our system consists of a chain of atoms interacting with each other and with a scalar "electromagnetic" field. 
Atoms are represented by their dipole moments P n , where the index n represents the position of an atom in the chain. 
Dynamics of the atoms is described within the tight-binding approximation with an interaction between nearest 
neighbors only, 

{Q 2 n - u 2 )P n + $(P n+1 + P„_i) = aE(x n ), (1) 

where $ is a parameter of the interaction, and f2 2 represents the site energy. Impurities in the model differ from host 
atoms in this parameter only, so 

n^ = ng Cn + n?(i-c l ) ) (2) 

where fig is the site energy of a host atom, flf describes an impurity, c„ is a random variable taking values 1 and 
with probabilities 1 — p and p, respectively. Parameter p, therefore, sets the concentration of the impurities in our 
system. This choice of the dynamical equation corresponds to exciton-like polarization waves. Phonon-like waves can 
be presented in a form that is similar to Eq. (|l|) with fl^ = Qq + (1 — c„)(l — Aide/ /Mhost)^ 2 , where M c i e f and Mh os t 
are masses of defects and host atoms, respectively. 

Polaritons in the system arise as collective excitations of dipoles (polarization waves) coupled to the electromagnetic 
wave, E(x n ), by means of a coupling parameter a. The electromagnetic subsystem is described by the following 
equation of motion 



c 2 ax' c 



where the right hand side is the polarization density caused by atomic dipole moments, and c is the speed of light in 
vacuum. The coordinate x in Eq. (|^) is along the chain with an interatomic distance a. Eqs. ([!]) and (^|) present a 
microscopic description of the transverse electromagnetic waves propagating along the chain in the sense that it does 
not make use of the concept of the dielectric permeability, and takes into account all modes of the field including 
those with wave numbers outside of the first Brillouin band. 

This approach enables us to address several general questions. A local state is usually composed of states with 
all possible values of wave number k. States with large k cannot be considered within a macroscopic dielectric 
function theory, and attempts to do so lead to divergent integrals that need to be renormalized.113 In our approach, 
all expressions are well defined, so we can check whether a contribution from large k is important, and if the long 
wave approximation gives reliable results. Calculation of the integrals appearing in the 3-1? theory requires detailed 
knowledge of the spectrum of excitations of a crystal throughout the entire Brillouin band. This makes analytical 
consideration practically unfeasible. In our 1-D model, we can carry out the calculations analytically (in a single- 
impurity case) and examine the influence of different factors (and approximations) upon the frequency of a local state 
and the transmission coefficient. Using caution, the results obtained can be used to assess approximations in 3-D 
cases. 
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2. A single impurity problem 



The equation for the frequency of the local polariton state in the 1-D chain has the form similar to that derived in 
Ref. [ i 



1 = Afi 2 G(0), 

where, however, the expression for the polariton Green's function G(n - 
of the system can be obtained in the explicit form 



(4) 



no) responsible for the mechanical excitation 



G(n - n ) = ■ 



cos(afc) - cos(^) 



- fig - 2$ cos (ka)] [cos(afc) - cos(^)] 
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exp [ik (n — no) a] 



(5) 



If one neglects the term responsible for the coupling to the electromagnetic field, the Green's function G{n — no) is 
reduced to that of the pure atomic system. This fact reflects the nature of the defect in our model: it only disturbs 
the mechanical (not related to the interaction with the field) properties of the system. A solution of Eq. (Q) can 
be real- valued only if it falls into the gap between the upper and lower polariton branches. This gap exists if the 
parameter $ in the dispersion equation of the polariton wave is positive, and the effective mass of the excitations in 
the long wave limit is, therefore, negative. 

The diagonal element, G(0), of Green's function (||) can be calculated exactly. The dispersion equation ||) than 
takes the following form: 



Afi 2 
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Qi,a(w) 



D(u>) 



.aw. to 2 — fin" 1 



± 2^H, 



,dUI . U> 2 — fin 

cos ( — ) nx 

v c ' 2$ 



A-KOLOJ CLLU 

— j— sm( — 

•Pc c 



give the poles of the integrand in Eq. (^|). The bottom of the polariton gap is determined by the condition D{u>) 
yielding in the long wave limit, auj/c <C 1, for the corresponding frequency, u>i, 



n z - 2%d 



(6) 

(7) 

(8) 
= 0, 

(9) 



where we introduce parameters d 2 = Ana/a, fip = fi§ + 2$, and take into account that the band width of the 
polarization waves, <£>, obeys the inequality y/^a/c <C 1. The last term in this expression is the correction to the 
bottom of the polariton gap due to the interaction with the transverse electromagnetic field. Usually this correction is 
small, but it has an important theoretical, and, in the case of strong enough spatial dispersion and oscillator strength, 
practical significance .El Because of this correction the polariton gap starts at a frequency, when the determinant D(u>) 
becomes imaginary, but functions Qi^tw) are still less than 1. This leads to the divergence of the right-hand side of 
Eq. (||) as to approaches a;;, and, hence, to the absence of a threshold for the solution of this equation. This divergence 
is not a 1-D effect since the same behavior is also found in 3-D isotropic modelflliJ An asymptotic form for Eq. (||) 
when uj — ► u>i in the 1-D case reads 



Afi 2 



and differs from the 3-D case by the factor of (diuia) / (cy/^j . The upper boundary of the gap, 
the condition Q\(lo) = 0, leading to 



'up ; 



(10) 



is determined by 



fin. 



(11) 
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Eq. (^J) also has a singularity as ui — > ui up , but this singularity is exclusively caused by the 1-D nature of the system. 
We will discuss local states that are not too close to the upper boundary in order to avoid manifestations of purely 
1-D effects. 

For frequencies deeper inside the gap, Eq. (^) can be simplified in the approximation of small spatial dispersion, 
y/Wa/c < 1, to yield 



u> 



ftf - Aft 2 



o 2 

"0 



ft 2 + 4$ 



Aft 2 



2c 



2 - ft 2 (ft 2 +(P -LO 2 



(12) 



where ft 2 = ft? + 2$ is a fundamental (k — 0) frequency of a chain composed of impurity atoms only. Two other 
terms in Eq. ( |12| ) present corrections to this frequency due to the spatial dispersion and the interaction with the 
electromagnetic field respectively. One can see that both corrections have the same sign and shift the local frequency 
into the region between ft 2 , and ftf . As we see below, this fact is significant for the transport properties of the chain. 

Transmission through the system can be considered in the framework of the transfer matrix approach. This method 
was adapted for the particular case of the system under consideration in Ref. ^. The state of the system is described 
by the vector v n , with components P n , P n +i, E n , E^/k^), which obeys to the following difference equation 



Un+l 



(13) 



The transfer matrix T n describes the propagation of the vector between adjacent sites: 
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T n = 



-1 -- 
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cos ka sin ka 

— sin ka cos ka / 



(14) 



Analytical calculation of the transmission coefficient in the situation considered is not feasible even in the case of a 
single impurity because the algebra is too cumbersome. The problem, however, can be simplified considerably if one 
neglects the spatial dispersion of the polarization waves. In this case the T-matrix can be reduced to a 2 x 2 matrix 
of the following form 



cos ka sin ka 

sin ka + f3 n cos ka cos ka + (3 n sin ka 



where the parameter (3 n 



fin 



c(c 2 -ft 2 )' 



(15) 



(16) 



represents the polarizability of the n-th atom due to its vibrational motion. In this case the complex transmission 
coefficient, t, can be easily expressed in terms of the elements of the resulting transfer matrix, T^ N ^ = ]X r„, 



ikL 



(rw + rw)-i(: 



T, 



(AO 



T. 



(N) 



(17) 



The problem is, therefore, reduced to the calculation of T^ N \ In the case of a single impurity, the product of the 
transfer matrices, r, can be presented in the following form 



1 K 1 = T "X T de f X 



-no — 1 



(18) 



where the matrix Tdef describes the impurity atom with ft„ = fti. The matrix product in Eq. ( jig ) is conveniently 
calculated in the basis, where the matrix t is diagonal. After some cumbersome algebra, one obtains for the complex 
transmission coefficient: 



t = 



2e lkL exp(-KL) 



1 - -fy (2 - P cot ka) [(l + e)]+ 2i exp (-kL)T cosh [na(N - 2n + 1] 



(19) 
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where R = /? 2 +4/?cot(a/c)-4, T = e(3/(sin{ka)VR), k is the imaginary wave number of the evanescent electromagnetic 
excitations, which determines the inverse localization length of the local state, and e = (fidef — P) /2yR. The last 
parameter describes the difference between host atoms and the impurity, and is equal to 

- (n?-ng) 



r (^-^K^-fif)' 

We have also neglected here a contribution from the second eigenvalue of the transfer matrix, which is proportional 
to exp(— 2kL), and is exponentially small for sufficiently long chains. For e = 0, Eq. ( |Tgj ) gives the transmission 
coefficient, to, of the pure system, 

2e ifeL e xp(-«L) 

1- (2-/3 cot fca)' 1 J 

exhibiting a regular exponential decay. At the lower boundary of the polariton gap, Oo, parameters [3 and k diverge, 
leading to vanishing transmission at the gap edge regardless the length of the chain, ft is instructive to rewrite Eq. 
( |l9| ) in terms of to'- 

*° (22) 



(1 + e) + i exp (-ikL) Tt cosh [na(N - 2n + 1] 

This expression describes the resonance tunneling of the electromagnetic waves through the chain with the defect. 
The resonance occurs when 

l + e = 0, (23) 
the transmission in this case becomes independent of the system size. Substituting the definition of the parameter 



e given by Eq. (|2C|) into Eq. (23), one arrives at an equation identical to Eq. (|12j) for the frequency of the local 
polariton state with the parameter of the spatial dispersion, being set to zero. The transmission takes a maximum 
value when the defect is placed in the middle of the chain, N — 2uq + 1 = 0, and in this case 

|W| 2 = ^ < 1. (24) 
The width of the resonance is proportional to Tto and decreases exponentially with an increase of the system's size. 



In the long wave limit, ai « 1, Eq. (24) can be rewritten in the following form 



(25) 



where u> r is the resonance frequency satisfying Eq. (|23j). ft is interesting to note that the transmission coefficient 
becomes exactly equal to one if the resonance frequency happens to occur exactly in the center of the polariton gap. 
This fact has a simple physical meaning. For uj 2 = fig + d 2 /2 the inverse localization length k becomes equal to the 
wave number uj r /c of the incoming radiation. Owing to this fact, the field and its derivative inside the chain exactly 
match the field and the derivative of the incoming field as though the optical properties of the chain are identical to 
those in vacuum. Consequently, the field propagates through the chain without reflection. 

Having solved the transmission problem we can find the magnitude of the field inside the chain in terms of the 
incident amplitude E^ n at the resonance frequency. Spatial distribution of the field in the local polariton state can 
be found to have the form E — Ed exp (— |n — n^na). Matching this expression with the outcoming field equal to 
Ei n texp(ikL) one has for the field amplitude at the defect atom, Ed, 

Ed = E in t exp(—ikL) exp [(N — no)na] . (26) 

For \t\ being of the order of one in the resonance this expression describes the drastic exponential enhancement of the 
incident amplitude at the defect side due to the effect of the resonance tunneling. 



Equations (22) and (|2J) demonstrate that the resonance tunneling via local polariton states is remarkably different 
from other types of resonance tunneling phenomena, such as electron tunneling via an impurity state,ll3 or through 
a double barrier. The most important fact is that the frequency profile of the resonance does not have the typical 
symmetric Laurentian shape. At uj — fij the parameter e diverges causing the transmission to vanish. At the same 
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time the resonance frequency u> r is very close to fix as it follows from Eq. (|12|). This results in strongly asymmetric 
frequency dependence of the transmission, which is skewed toward lower frequencies. 

The transmission vanishes precisely at two frequencies: at the low frequency band edge Jlo and at the fre- 
quency f2i associated with the vibrational motion of the defect atom. At the same time, the behavior of the 
transmission coefficient in the vicinities of these two frequencies is essentially different: at the band edge it is 

(lu 2 — f2Q) 2 cxp ^— 1/ yjuj 2 — £1q^ while at the defect frequency the transmission goes to zero as (lu 2 — ^ 2 ) 2 - These 

facts can be used to predict several effects that would occur with the increase of the concentration of the defects. 
First, one can note that with the increase of concentration of the impurities frequency fl\ becomes eventually the 
boundary of the new polariton gap when all the original host atoms will be replaced by the defects atoms. One can 
conclude then that the zero of the transmission at fii instead of being washed out by the disorder, would actually 
become more singular. More exactly one should expect that the frequency dependence of the transmission in the 
vicinity of f2i will exhibit a crossover from the simple power decrease to the behavior with exponential singularity 
associated with the band edge. Second, if one takes into account such factors as spatial dispersion or damping, which 
prevent transmission from exact vanishing, one should expect that the above mentioned crossover to the more singular 
behavior would manifest itself in the form of substantial decrease of the transmission in the vicinity of f^i with an 
increase of the concentration. Numerical calculations discussed in the next section of the paper show that this effect 
does take place even at rather small concentration of the defects. 

Resonance tunneling is very sensitive to the presence of relaxation, which phenomenologically can be accounted for 
by adding 2i"/Lu to the denominator of the polarizability /3, where 7 is an effective relaxation parameter. This will 
make the parameter e complex valued, leading to two important consequences. First, the resonance condition becomes 
Re(e) = — 1, and it can be fulfilled only if the relaxation is small enough. Second, the imaginary part of e will prevent 
the exponential factor to in Eq. (^2|) from canceling out at the resonance. This restricts the length of the system in 
which the resonance can occur and limit the enhancement of the field at the defect. These restrictions though are not 
specific for the system under consideration and affect experimental manifestation of any type of resonant tunneling 
phenomenon. 

Since we only concern with a frequency region in the vicinity of Q±, real, £1, and imaginary, £2, parts of £ can be 
approximately found as 

2 n ia I aTj 2 u 2 -n 2 

^\ d 2 -An 2 + {7 > 



£ ^Z^k £l - (28) 

It follows from Eq. ( [27j ) that the resonance occurs only if (47c) /(ad 2 ) < 1. This inequality has a simple physical 
meaning: it ensures that the distance between the resonance frequency, u> r , and Qi, where the transmission goes 
to zero, is greater than the relaxation parameter, 7. This is a rather strict condition that can only be satisfied for 
high frequency oscillations with large oscillator strength in crystals with large interatomic spacing, a. The spatial 
dispersion, however, makes conditions for the resonant tunneling much less restrictive. In order to estimate the effect 
of the dissipation in the presence of the spatial dispersion one can rely upon Eq. ( |22] ) assuming that the dispersion 
only modifies the parameter s, but does not effect the general expression for the transmission. This assumption is 
justified by the numerical results of Ref. ^ and the present paper, which show that the transmission properties in 
the presence of the spatial dispersion do not differ significantly from the analytical calculations performed for the 
chain of noninteracting dipoles. According to Eq. (p"2|), the inter-atomic interaction moves the resonance frequency 
further away from fij undermining the influence of the damping and leading to a weaker inequality: (7^1)/$ < 1. 
This condition can be easily fulfilled, even for phonons with a relatively small negative spatial dispersion. For the 
imaginary part £2 at the resonance one can obtain from Eq. (Eq) the following estimate 



£2 ~ min[(47c)/H 2 ), ( 7 fii)/*]. (29) 

The requirement that £2 be much smaller than to leads to the following restriction for the length of the system 
L <C (1/k) I ln[£2] I, with £2 given above. The maximum value of the field at the defect site attainable for the defect 
located in the center of the chain is then found as \Ed\ ~ \Ei n \\t\/ y/e2. 
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III. ONE-DIMENSIONAL DIPOLE CHAIN WITH FINITE CONCENTRATION OF IMPURITIES 



In this section we present results of numerical Monte-Carlo simulations of the transport properties of the system 
under consideration in the case of randomly distributed identical defects. If spatial dispersion is taken into account 
the regular Maxwell boundary conditions must be complemented by additional boundary conditions regulating the 
behavior of polarization P at the ends of the chain. In our previous papercl we calculated the transmission for two types 
of boundary conditions: Pq = Pn = 0, which corresponds to the fixed ends of the chain, and Pq = Pi, fV-i = Pn, 
which corresponds to the relaxed ends. We reported in Ref. ^ that the transmission is very sensitive to the boundary 
conditions with fixed ends being much more favorable for the resonance. Our present numerical results obtained 
with an improved numerical procedure and the analytical calculations do not confirm this dependence of the resonant 
tunneling upon the boundary conditions. In the case of a single impurity we find that for both types of the boundary 
conditions the transmission demonstrates sharp resonance similar to that found in Ref. ^ for fixed ends. Similarly, for 
a finite concentration of impurities we did not find any considerable differences in the transmission for both types of 
boundary conditions. We conclude that the actual form of the boundary conditions is not significant for the resonant 
tunneling. 

The transfer matrix, Eq. (|l3|), along with the definition of the transfer matrix, Eq. (p"4"|), and the boundary 
conditions chosen in the form of fixed terminal points, provides a basis for our computations. However, it turns out 
that straightforward use of Eq. (|l3|) in the gap region is not possible because of underflow errors arising when one 
pair of eigenvalues of the transfer matrix becomes exponentially greater than the second one. In order to overcome 
this problem we develop a new computational approach based upon the blend of the transfer-matrix method with the 
invariant embedding ideas. The central element of the method is a 4 x 4 matrix S(N) that depends upon the system 
size N. The complex transmission coefficient, t, is expressed in terms of the elements of this matrix as 

t = 2 ex-p(-ikL)(S n + S 12 ). (30) 

The matrix S(N) is determined by the following nonlinear recursion 

S(N + 1) = T N x S(jV) x S(N), (31) 

where matrix 3(iV) is given by 

E(N) = {I-S(N) x H x [I-T{N)}}-\ (32) 
The initial condition to Eq. (|3l]) is given by 

5(0) = (G + H)- 1 , (33) 

where matrices G and H are specified by the boundary conditions. The derivation of the Eqs. (|30|)-(|33l) and more 
detailed discussion of the method is given in the Appendix to the paper. The test of the algorithm based upon 
recursion formula (|3l| ) proves the method to provide accurate results for transmission coefficients as small as I0 -15 . 

In our simulations we fix the concentration of the defects and randomly distribute them among the host atoms. 
The total number of atoms in the chain is also fixed; the results presented below are obtained for a chain consisting of 
1000 atoms. For the chosen defect frequency, Sli ~ 1.354f2 , the localization length of the local polariton state, k n d,is 
approximately equal to 150 interatomic distances. The transmission coefficient is found to be extremely sensitive to 
a particular arrangements of defects in a realization exhibiting strong fluctuations from one realization to another. 
Therefore, in order to reveal the general features of the transmission independent of particular positions of defects, we 
average the transmission over 1000 different realizations. We have also calculated the averaged Lyapunov exponent 
(the inverse localization length, I chain, characterizing transport through the entire chain) to verify that the averaged 
transmission reveals a reliable information about the transport properties of the system. 

The results of the computations are presented in the figures below. Figs. 1-3 show an evolution of the transmission 
with the increase of the concentration of the impurities. In Fig. 1 one can see the change of the transport properties 
at small concentrations up to 1%. The curve labeled (1) shows, basically, the single impurity behavior averaged over 
random positions of the defect. With an increase of the concentration there is a greater probability for two (or more) 
defects to form a cluster resulting in splitting a single resonance frequency in two or more frequencies. The double 
peak structure of the curves (2) and (3) reflects these cluster effects. With the further increase of the concentration the 
clusters' sizes grow on average leading to multiple resonances with distances between adjacent resonance frequencies 
being too small to be distinguished. Curve (5) in Fig. 1 reflects this transformation, which marks a transition between 
individual tunneling resonances and the defect induced band. The concentrations in this transition region is such 
that an average distance between the defects is equal to the localization length of the individual local states, k n d- 
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The collective localization length at the frequency of the transmission peak, l^^ in , becomes equal to the length of the 
chain at approximately the same concentration that allows to suggest a simple relationships between the two lengths: 
I chain = °hnd, where c stands for the concentration. The numerical results presented in Fig. 4 clearly demonstrate this 
linear concentration dependence of l™j^ in at small concentrations. For larger concentrations one can see from Figs. 2 
and 3 that the peak of the transmission coefficient developes into a broad structure. This marks further development 
of the defect pass band. Curves in Fig. 2 show the transmission coefficient at intermediate concentrations, where 
localization length, l c hain, is bigger than the length of the system only in a small frequency region around the maximum 
of the transmission, and Fig. 3 presents well developed pass band with multipeak structure resulting from geometrical 
resonances at the boundaries of the system. 

This figures reveal an important feature of the defect polariton band: its right edge does not move with increase of 
the concentration. The frequency of this boundary is exactly equal to the defect frequency, (which is normalized 
by Qq in the figures), and the entire band is developing to the left of S~2i in complete agreement with the arguments 
based upon analytical solution of the single-impurity problem. Moreover, the magnitude of the transmission in the 
vicinity of ^decreases with increase of the concentration also in agreement with our remarks at the end of the previous 
section. Fig. 5 presents the inverse localization length, l c hain, normalized by the length of the chain for three different 
concentrations. It can be seen that l~u ain { significantly grows with increase of the concentration, reaching the 
value of approximately 17 /L at the concentration as small as 3%. Such a small localization length corresponds to the 
transmission of the order of magnitude of 10 -17 , which is practically zero in our computation. Further increase of the 
concentration does not change the minimum localization length. These results present an interesting example of the 
defects building up a boundary of the new forbidden gap. 

From this figure it is also seen the development of the pass band to the left of Vt\ presented above in Figs. 1-3, but 
at a larger scale. We can not distinguish here the details of the frequency dependence, but the transition from the 
single resonance behavior to the pass band, marked by the significant flattening of the curve, is clear. 

The last Fig. 6 presents the concentration dependence of the semiwidth, 6 U , of the defect band. The semiwidth 
is defined as the difference between frequency of the maximum transmission and the right edge of the band. One 
can see that all the point form a smooth line with no indication of a change of the dependence with the transition 
between different transport regimes. Attempts to fit this curve showed that it is excellently fitted by the power law 
S u cx c^with v ~ 0.8 in all studied concentration range. The reason for this behavior and why it is insensitive to the 
change of the character of the transport requires further study. 

IV. CONCLUSION 

In this paper, we considered one-dimensional resonance tunneling of scalar "electromagnetic waves" through an 
optical barrier caused by a polariton gap. The tunneling is mediated by local polariton states arising due to defect 
atoms embedded in an otherwise ideal periodic chain. We also numerically studied how a defect-induced propagating 
band emerges from these resonances when the concentration of defects increases. 

It is important to emphasize the difference between the situation considered in our paper and other types of 
tunneling phenomena discussed in the literature. The tunneling of electromagnetic waves through photonic crystals 
and electron tunneling, despite all the difference between these phenomena, share one common feature. In both cases, 
the resonance occurs due to defects that have dimensions comparable with wavelengths of the respective excitations 
(electrons interact with atomic impurities, and long wave electromagnetic waves interact with macroscopic distortions 
of the photonic crystals). In our case the wavelength of the propagating excitations is many orders of magnitude 
greater than dimensions of the atomic defects responsible for the resonance. The physical reason for such an unusual 
behavior lies in the nature of local polaritons. These states are formed owing to the presence of internal polariton- 
forming excitations. The spatial extent of these states is much larger than the geometrical dimensions of atomic 
defects and is comparable to the wavelength of the incident radiation. 

We presented an exact analytical solution of the tunneling of electromagnetic waves through a chain of noninteracting 
atoms with a single defect. This solution provides insight into the nature of the phenomenon under considerarion 
and allows one to obtain an explicit expression for the magnitude of the electromagnetic field at the defect site. The 
expression derived demonstrates that the field is strongly enhanced at the resonance with its magnitude growing 
exponentially with an increase of the length of the system. This effect is an electromagnetic analogue of the-,diarge 
accumulation in the case of electron tunneling, where it is known to cause interesting nonlinear phenomena.li3c3 

An analytical solution of the single-defect problem allowed us to make predictions regarding the transport properties 
of the system with multiple randomly located defects. The most interesting of these is that the dynamical frequency 
of the defects, fix, sets a high frequency boundary for the defect induced pass band, which does not move with 
increasing concentration of defects. Numerical Monte-Carlo simulations confirmed this assumption and showed that 
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the direct interaction between the atoms (spatial dispersion) does not affect resonance tunneling considerably, though 
it adds new interesting features to it. One of them is the behavior of the transmission in the vicinity of Cl±. In 
absence of the spatial dispersion, the transmission at this point is exactly equal to zero, and remains small when the 
interaction is taken into account. The interesting fact revealed by the numerical analysis is that the transmission at 
Q\ decreases with an increase in the concentration of the defects and nearly approaches zero at concentrations as small 
as 3%. This fact can be understood in light of the transfer matrix approach: if the frequency f2i corresponds to the 
eigenvalue of the defect's transfer matrix, which significantly differs from one, the transmission will diminish strongly 
each time the wave encounters a defect site, regardless the order in which the defects are located. Numerical results 
also demonstrated a transition between two transport regimes: one associated with resonance tunneling and the other 
occurring when the resonances spatially overlap and a pass band of extended states emerges. The transition occurs 
when the average distance between the defects becomes equal to the localization length of the single local state. At the 
same time the collective localization length at the peak transmission frequency, characterizing the transport properties 
of the entire chain, becomes equal to the total length of the system. This result assumes the linear dependence of 
this collective localization length upon concentration, which we directly confirm for small concentrations. Numerical 
results also showed that the width of the resonance, which develops into a pass band with an increase in concentration, 
does not manifest any transformation when the character of transport changes. The concentration dependence of the 
width was found to be extremely well described by a power law with an exponent approximately equal to 0.8. The 
nature of this behavior awaits an explanation. 



APPENDIX A: INVARIANT EMBEDDING ALGORITHM FOR THE TRANSFER-MATRIX EQUATION 



In this Appendix we develop invariant embedding approach to transfer-matrix equations of a general form and 
deduce Eqs. (|3^)-(|3^) used for Monte-Carlo calculations in our paper. Consider a typical difference equation of the 
transfer-matrix method, 



with boundary conditions of a general form: 



^n+l — T n ?i n , 



Guq + Hun = v. 



(Al) 



(A2) 



Here u n is a vector of an appropriate dimension that characterizes the state of the system at the nth site, T n is a 
respective transfer-matrix; G and H are matrices of the same dimension as the transfer-matrix, together a with the 
vector v they specify boundary conditions at the left and right boundaries of the system (cites and N respectively). 
The regular Max wel l boundary conditions and the fixed ends boundary condition for polarization can be presented 
in the form Eq. (A2) with the following matrices G,H, and vector v 



G 



1 


—i 








1 


—i 














1 











1 






H = 



1 i 

-1 -i 

1 

-1 



(A3) 



These matrices are singular, but one should not worry about this, because we will only neeeLto invert their sum, which 
has a non-zero determinant. In accord with the ideas of the invariant embedding methodll3 we consider the dynamic 
vector, u n , as a function of the current site n, the length of the system N, and the boundary vector v: 



u(n,N,v) = S(n,N)v 



(A4) 



In the last equation we use the line ar n ature of E q. (Al) in order to separate out the dependence upon the vector 
v. Substituting Eq. (AA) into Eqs. (|Al[) and (A2) we have the dynamical equation and boundary conditions for the 
matrix S : 



S(n + 1, N) = T n x S(n, N), 



(A5) 



G x 5(0, N) + H x S(N, N) = I, 



(A6) 



where / is a unit matrix. The matrix S(n, iV + 1), which describes the system with one additional scatterer, obviously 
satisfies the same equation ( |A5| ) as S(n,N). Relying again upon the linearity of Eq. (A5) we conclude that S(n,N) 
and S(n,N + 1) can only differ by a constant (independent of n) matrix factor A(N). 
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S(n, N+l) = S{n, N) x A(N). 



(A7) 



In order to find A(N) we first substitute Eq. (A7) into boundary conditions Eq. (A6) which yield 

A(N) = Gx S(0, N + 1) + H x S(N, N + 1). 
Boundary conditions Eq. (|A6|) do not change if A is replaced by A + 1 , therefore we can write down that 

Gx 5(0, A + 1) = I — H x S(N + 1,N+1). 



Substituting this expression into Eq. ( A8) we have for the matrix A(N) : 

A(N) = I + H x [S(N, N + 1) - S(N + 1,N+ 1)]. 



(A8) 



(A9) 



(A10) 



The quantity S(A + 1, TV + 1) can be eliminated from this equation by means of Eq. (Al): S(N + 1, N + 1) = Tjy 
S(N, N + 1), and we have for A(A) 



A(N) =I + Hx[I- T{N)} x S(A, A + 1). 



(All) 



Substituting this formula into Eq. ( A7) we obtain the equation that governs the evolution of the matrix S(n, N) with 
the change of the parameter A: 



S(n, N + 1) = S(n, A) + S(n, N) x H x [I - T(N)} x S(N, N + 1). 



(A12) 



This equation, howe ver, is not closed because of an unknown matrix S(N, N + 1). This matrix can be found by setting 
ri = N in Eq. (|A12|) : 



S(N, N+1) = {I- S(N, N)xHx[I- r(JV)]} _1 5(iV, N). 

Introducing notation 

S(N) ={I- S(N, N)x H x[I - T(N)]}~ 1 
the previous expression can be rewritten in the following compact form: 

S(N, N+l) = S(N) x S{N, N). 
Inserting Eq. ( A15| ) into Eq. ( Al2j ) we finally obtain: 

S(n, N + l) = S(n, N) + S(n, N) x H x [I - T(N)] x S(N) x S(N, N). 



(A13) 
(A14) 
(A15) 
(A16) 



This equation still has an unknown qu antit y S(N, N) whic h must be determined separately. We achieve this combining 
the original transfer matrix equation ( |Al[ ) and Eq. (A15) to obtain the following: 

(A17) 



S(N + 1, N + 1) = T N x E(N) x S(N, N). 



Eq. (A17) is a nonlinear matrix equation with an initial condition given by 

(G + H)x 5(0,0) = L 



(A18) 



Eqs. (|30|)-(|33|) of the main body of the paper coincide with Eqs.( [A16| )-( A18 ) with simplifyed notation for the matrix 
S, where we dropped the second argument. They constitute the complete set of embedding equations for the transfer 
matrix problem. In order to find the transmission coefficient one has to multiply the matrix S(N, N) by the boundary 
vector v; the first component of the resulting vector is equal to texp(ikL), where t is the complex transmission 
coefficient. If one is interested in the d istrib ution of the state vector u(n, N) throughout the entire system, one has 
to find S(N,N) and then to solve Eq. (|Al6j) . 

The presented algorithm was proved to be extremely stable, it produced reliable results for tr ansm ission as small as 
10 -17 . This stability is due to the operation of inversion involved in the calculations [see Eq. (A14)]. This operation 
prevents elements of the matrix S to grow uncontrollably in the course of calculations. 
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Figure Captions 

Fig. 1. Frequency dependence of the averaged transmission coefficient for small concentrations of the defects. The 

frequency is normalized by the fundamental (k = 0) frequency of the pure chain, Slo- The low- frequency boundary of 

the polariton gap is at u » 1.3 and is not shown here. 

Fig. 2. The same as in Fig. 1 but for intermediate concentrations. 

Fig. 3. The same as in Fig. 1 but for large concentrations. 

Fig. 4. Concentration dependence of the collective localization length, l c hain, normalized by the system's size L. 
Fig. 5. Frequency dependence of the Lyapunov exponent of the entire chain for several concentrations in the frequency 
region of the defect band. 

Fig. 6. Concentration dependence of the semiwidth of the defect band. The solid line represents fit with power 
function c v ', where v ss 0.8. 
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